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Abstract 



^ ■ The discrete Gaussian model for the surface of a crystal deposited on a 

jrt ' disordered substrate is studied by Monte Carlo simulations. A continuous 

transition is found from a phase with a thermally-induced roughness to a 

C . glassy one in which the roughness is driven by the disorder. The behavior of 

CJ ■ 

the height-height correlations is consistent with the one-step replica symmetry 

^ . broken solution of the variational approximation. The results differ from the 

renormalization group predictions and from recent simulations of a 2D vortex- 
glass model which belongs to the same universality class. 
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The surface of a crystal formed by the deposition of particles on a flat substrate is 
undergoing a roughening transition: At temperatures below T^, the roughening transition 
temperature, the surface is flat. For T > Tr the periodic potential due to the discrete and 
uniform size of the particles is irrelevant and the thermal fluctuations cause the surface to 
roughen [Q,0. The intriguing question of what will be the crystalline surface properties 
with a disordered substrate has been addressed recently |^. The renormalization group 
(RG) analysis predicts a low-temperature superrough phase with surprisingly uncommon 
equilibrium and dynamic |P3 characteristics. Near the superroughening transition the 
Hamiltonian of the surface is that of the 2D sine-Gordon model with random phases |^|-0. 
Another system which belongs to the same universality class (and therefore should follow a 
similar behavior near criticality) is the two dimensional array of flux lines with the magnetic 
field parallel to the superconducting plane in the presence of random pinning ||^ . In that 
context the random sine-Gordon model (RSGM) has been studied by a non-perturbative 
variational (harmonic) approximation [p|-[Tl|] which allows replica-symmetry breaking with 
a low-temperature phase drastically different than that found by the RG approach |]T^,|I3 



Recent numerical simulations in the weak coupling limit of the RSGM [|T1| have not observed 
any signature of a transition in the equilibrium properties. A transition was found, however, 
in the dynamic response in qualitative agreement with the dynamic RG |^. 

In view of these results we present here the study by large scale Monte Carlo (MC) 
simulations of the disordered-substrate surface model. The main conclusion is the exis- 
tence of a low-temperature glassy rough phase separated by a continuous transition (at 
the critical temperature anticipated by the analytic calculations) from the thermally-rough 
high-temperature phase. Height-height correlations in the glassy phase are consistent with 
the predictions of the one-step replica symmetry broken solution obtained from the varia- 
tional approach [jT2|,|13[. Finite-size scaling for the ratio of moments of a "modified order- 
parameter" provides a first rough estimate for the correlation length exponent at the tran- 
sition. 

The surface configurations are described in terms of height variables hi, defined on every 



lattice point i of the 2D basal plane. The Hamiltonian we choose is that of the discrete 
Gaussian model: 
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where k is the surface tension. In the case of a fiat surface the hi takes integer values in 
the units of the lattice spacing a in the direction perpendicular to the surface. To simulate 
the disordered substrate we first assign a random quenched height di, chosen uniformly (and 
independently for each site), in the interval (— a/2, +a/2]. The height hi is then forced to 
take the values hi = di + riia where rii is any, positive or negative, integer. 

The calculation of the partition function for a given realization of the disorder consists 
in summing exp {— /37i} over all integers rij. The relation to the RSGM is obtained by using 
the Poisson summation formula. Introducing continuous fields 0j, the partition function 
may be expressed as: 
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In the continuum limit 0j are then replaced by 0(x) and the discrete Laplacian by a contin- 
uous one. In addition, near the critical point only the first harmonic in the local periodic 
"potential" for cj) is relevant ||I|. When higher (irrelevant) harmonics are discarded the 
Hamiltonian of the RSGM is obtained: 



Z = j d(f){x) exp <— dx 



— [V0(f)]' - A cos (27r [0(f) - d{x)] /a) 



(3) 



The equilibrium scaling properties are manifested in the behavior of the height-height 
correlation function: 



C{v)= Uh{v + vo)-h{vo)y 



(4) 



where (■ ■ ■)t denotes a thermal (time) average for a given realization of disorder and [■ ■ ■]av 
represents the configurational average over different realizations of the disorder. Overbar 
denotes the averages over all origins of tq and all directions. 



We briefly recall the major results previously obtained for the C(r): Analytically the 
Renormalization group |p and the variational [|13,|1^ calculations both predict a transition 
at the same critical temperature Tc = n/n. Both approaches also agree on the properties of 
the high temperature phase in which the discrete nature of the particle is irrelevant and the 
behavior is that of the simple Gaussian model [without the cosine term in the Hamiltonian 
of Eq. d^)]. The behavior of C(r) is therefore given by: 

C(r) = — ln|r|. (5) 

nil 

As mentioned above the two approaches diverge in their predictions regarding the be- 
havior for T < Tc. RG calculations PJ^,^,0 predict a new term with (ln|r|)^ which will 
dominate at large distances. More precisely the first-order calculations yield the following 
behavior: 

Cir) = ACoir) + BT\\n\r\f, (6) 

where r = 1 — T/Tc, Co(r) = T/{ttk,) In |r| is the correlation function (^, A is a nonuniversal 
constant, and B = 2/tt'^ + 0{t) is a universal constant. 

The variational approach [p!2| , p!3| , on the other hand, predicts the same scaling in In |r| 



as for T > Tc. The prefactor, however, becomes T-independent and sticks to its value at Tc 
for the whole T < Tc phase: 

C(r) = — Inlrl. (7) 

The MC calculations were performed using the Metropolis algorithm for the discrete 
Gaussian Hamiltonian Eq. (|ID on a square lattice with periodic boundary conditions. The 
simulations were carried on the CM computers of the Thinking Machines Corporation. At 
every time step all the variables hi of one sublattice were simultaneously updated by in- 
creasing or decreasing them (independently) by one unit. The moves are then accepted or 
rejected according to the Metropolis rules. 

Following the approach introduced by Young in spin-glass MC |T5|JT^, for every real- 
ization of disorder two replicas of the system were simulated. They had random initial 



conditions (realized by adding random integer heights on top of the random substrate) and 
each had its own independent time evolution. In addition to the information extracted from 
their overlap (see below) the two-replica approach allows close monitoring of the approach 
to equilibrium. After equilibration was established, the measurements were taken over a 
time interval which was one or several equilibration times long. Depending on the lattice 
size and the length of the MC runs, the disorder averaging were performed using 100 to 
several thousands realizations of the disorder. 

The Monte-Carlo simulations of the static height-height correlation function of the pure 
discrete Gaussian model were performed by W. Shugard et al [|l^. They confirmed the 
existence of the phase transition between a high-temperature rough phase characterized by 
a logarithmic behavior of C(r) and a fiat low-temperature phase. We reproduced their 
results as a special case of our model without disorder. However, introducing the disorder 
drastically changes the character of the low-temperature phase. 

Fig. |I| shows on a semilog plot the behavior of C(r) around the critical temperature 
Tc = k/tt = 0.6366 for our choice of k = 2. The simulations were performed with maximum 
lattice size of L = 64. For each temperature, measurements were repeated, with a new set of 
100 realizations of the disorder, five to ten times. The average values of these measurements 
with corresponding error bars are shown in Fig. ^ 

We compare our results with both the renormalization group predictions Eq. (P) and 
with the results derived by the variational analysis with one-step symmetry breaking scheme 
Eq. (^. In both cases the theoretical results refer to the large |r| behavior of C(r) while 
numerical results are always limited by the lattice size L and by finite-size effects. 

According to the RG (|^), for T near Tc and for large |r| the effect of the second term, 
i?r^(ln |r|)^, should dominate. To compare with Eq. (^) we need to be close to T^ but not 
too close (since the coefficient is proportional to r^ and the distance |r| at which this term 
will dominate will be beyond the size of the system). We therefore choose T = 0.5 for the 
comparison. The broken line shows Eq. (^ for that temperature and A = 1. It can be seen 
clearly that the upper bending trend of the broken line is inconsistent with the MC data. 



The down bending of the MC points for every temperature (which was also observed in the 
data taken for the pure system) are due to the finite-size effects. 

The rephca variational approach predicts (see Eq. (|^) C(r) to remain logarithmic for 
all T with a T-independent coefficient for T < T^. The data shown in Fig. |1| is consistent 
with this behavior. In our fit we used the values of the C(r) for |r| between 4 and 14 lattice 
spacings. We neglected the higher values of C(r) because of the presence of the finite size 
effects as well as of strong sample to sample fluctuations. The full straight lines in Fig. |l| 
are the best fit curves. In Fig. ^ we show the temperature dependence of the slopes for nine 
values of T between 0.45 and 0.9 including those in Fig. |I|. The vertical doted line T = k/tt 
is the analytic result for T^. While in high-T phase the slope of C(r) changes linearly with 
T, for the low-T phase it saturates around the value l/vr^ as is predicted by Eq. (|^). 

To gain more insight in the transition and the properties of the low-temperature phase, a 
glassy order-parameter, its correlations, and/or probability distribution, should be invoked. 
We tried to look at the local autocorrelation function: 

qfit) = { [htito + t)- h-{to + t)] [/if (to + eafst) - hl'ito + eafst)] } , (8) 

where to is an initial time (larger than the time required to equilibrate the system) and 
overbar means average over lattice sites. The replica indices a,f3 = take values 1 or 2 and 
ea/3 = if a = /5 or 1 if a 7^ /3. This local quantity is first averaged over all sites: 

1 ^ 
?"/3W = T7E?fW, (9) 

^^ i=i 

N = LxL. This order-parameter definition includes the equal time overlap between different 
replicas as well as the auto-overlap of the same replica at different times. In the limit of 
infinite time separation in the latter the probability distribution for both should coincide. 
It is Gaussian for T > T^ but is expected to deviate from it for T < T^. One can then look 
for the phase transition by calculating the ratio of the moments 



1 f. Ms) 



i3lT\av 



«'"'^>^5f-ii5i:l ''•" 



after thermal equilibrium has been reached. If the transition is of the second order the 
expected finite size scaling of this quantity is g{T,L) ~ g^L^/^iT — Tc)) where u is the 
correlation length exponent, ^ ^ \T — Tc\~'^ . 

Trying to evaluate qa/B and g we could not reach enough accuracy to extract reliable 
results within the computer time available. A similar problem has been observed in simula- 
tions of the 3D Ising spin glass [J16l and the 3D gauge glass ll8|. One possible approach is to 



look for a quantity with similar scaling but for which the statistical errors are smaller. The 
problem we have to overcome is that the deviations from the Gaussian distribution occur 
first (for T just below Tc) at very small values of q compared with its rms. To accentuate 
the contribution from these small values we tried to evaluate the "renormalized" quantity 

qap{t)- 



j:l,\qf{t)\ 



?-/3(^) = ^N I a/3.,M - (11) 



The distribution P{q) (see, e.g., Ref. |T6| for its definition) indeed exhibits a transition 
from a distribution with one maximum at g = to a distribution with two maxima symmetric 
with respect to g = (which becomes a local minimum). The figures of P{q) will be exhibited 



in a future publication |T^ . Here we only depict in Fig. ^ the result of the quantity g defined 



in Eq. (|10]) but with q replacing q. Equilibration was confirmed by the convergence of qaa 
and qai^. Depending on the temperature and the size of the simulated system, the MC runs 
for equilibration were between 2^^ to 2^^ MC steps. For smaller system size the average over 
disorder is performed using over 1000 samples. However, for large system size, L = 64, the 
equilibration time for temperatures below Tc is long, and we worked with a hundred samples. 
As in the calculation of C(r), for each temperature the measurements were repeated several 
times. Again, the average values of these measurements with corresponding error bars 
are shown in Fig. ^. Clearly, Fig. ^ strongly suggests existence of the phase transition 
somewhere in the temperature range 0.63-0.66. A finite size scaling plot of g is shown in 
inset of Fig. |. The best fit is obtained with Tc = 0.643 ± 0.006 and u = 1.23 ± 0.10. The 
value of Tc is in good agreement with the analytic predictions discussed above. 



We would like to emphasize that the main purpose in presenting the data for g is to 
provide an independent confirmation for the existence of the transition. We believe that 
replacing g by g will not change the critical temperature. In addition, usually the same 
exponent v controls the finite-size effects for all thermodynamic quantities, so the value 
found here may be a first rough estimate for this critical exponent. There is no analytic 
prediction for v within the variational approach [the RG approach [^] predicts ^ ~ exp(yl/r^) 
and is, again, in disagreement with our data]. Our plots of P{q) do not exhibit replica 
symmetry breaking but they are not accurate enough to rule it out. Simulations of larger 
systems should resolve this issue. 

Finally, we comment on the only other numerical work [Q, in which no transition was 



found in the equilibrium C(r): These simulations were performed in the weak coupling 
regime A <^ 1 of the DSGM. In contrast the ones presented here are in the strong coupling 
regime with all harmonics having coefficients 0(1). This raises the possibility that these 
cases belong to distinct universality classes. Whether such a scenario may be reconciled 
with the phase transition in the dynamics found in the same weak-coupling simulations [T3|, 
remains to be seen. 

To summarize, numerical evidence has been presented, from both the height-height cor- 
relations and the overlap distribution, for the existence of a phase transition in the scaling 
properties of a surface upon a disordered substrate. The low-temperature phase has glassy 
characteristics. The behavior of C(r) is consistent with that predicted from the one-step 
replica symmetry breaking solution although no direct evidence for (or against) such a 
breaking was found. Is this symmetry broken? Why there is such a discrepancy between the 
numerical results and these of the RG approach which was so successful in explaining the 
behavior of the pure (flat-substrate) surface? These are two among essential unanswered 
questions. More studies will be needed to reach a fuller understanding of the disordered- 
substrate surface, the fiux-lines array in 2D type-II superconductors, and other physical 
systems related to the random-phase sine-Gordon model. 
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FIGURES 
FIG. 1. Semilog plot of C(r) for L = 64. The straight lines are the best fits to 

C(r) = a{T) + b{T) In |r|. The broken line is the RG prediction Eq. (^) at T = 0.5 with A = 1. 

FIG. 2. Plot of the coefficient b{T) from the fitting equation C(r) = a{T) + b{T) ln|r|. The 
vertical dashed line is the analytic Tc. The horizontal line is the slope predicted by the Eq. (1^) for 
aU T < Tc. 

FIG. 3. Plot of g{T,L) vs temperature T for different lattice sizes. The critical temperature 
is indicated by the crossing of the curves. The full line is the guide to the eye. The inset shows 
the finite size scaling plot with the indicated values of Tc and v (see text). 
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